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ABSTRACT 


"Kalman filtering techniques utilizing a rotated coordi- 
nate system were applied to tracking problems encountered in 
fire control systems. Two models of target motion were con- 
sidered: a constant-velocity model and a model which assumes 
correlated random accelerations. Estimators derived from 
these models were evaluated using Monte-Carlo simulations of 
constant-velocity and maneuvering targets. An algorithm de- 
veloped to calculate prediction accuracy data for time in- 
tervals based on an approximation of the time of flight for 
a 5/54 Caliber projectile was used to obtain prediction 
accuracy statistics for evaluating estimator performance. 

A track generating program was developed to produce maneu- 
vering target state trajectories in either cartesian or 
spherical coordinates for use in the Monte-Carlo simulations. 
Adaptive filters were developed to take advantage of the 
best features of the two basic estimators and Monte-Carlo 


Simulation of the adaptive filters was performed. 
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I. INTRODUCTION 


A. ALTERNATIVE ESTIMATION TECHNIQUES 

Recently the Gun-Fire Control System MK 86 has been the 
subject of several studies done at the Naval Postgraduate 
School involving the use of Kalman filtering techniques in a 
gun-fire-control system [1],[2],[4]. This thesis examines 
some possible alternatives to the techniques used in GFCS 
MK 86 and makes comparisons to the GFCS MK 86 filter. 

One proposed alternative to GFCS MK 86 procedure is the 
use of a rotated cartesian coordinate system. The coordi- 
nate system is rotated about the vertical (z) axis so that 
the x axis is aligned with the initial measured target bear- 
ing. Figure 1 illustrates the proposed coordinate system 
and its relationship to the fixed cartesian and spherical 
coordinate systems when noiseless measurements are assumed. 
In the general case, when the measurements are noisy, the 
x-axis will not be aligned exactly with the true target 
bearing as indicated in Fig. 1 but the error is expected to 
be small. The advantages gained by rotating the coordinate 
system are diminished correlation of the measurement noise 
and more accurate modeling of random forcing inputs. These 
attributes are discussed quantitatively in Chapter II. 

In Chapter III estimators based on two models of target 
motion are developed: a constant-velocity model similar to 
one used іп GFCS MK 86 (1] but utilizing the rotated coordi- 


nate system, and a constant-velocity model with a coloring 
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pre-filter similar to that described by McKinley [4] also 
using the rotated coordinate system. 

Monte-Carlo simulations were performed in order to eval- 
uate the performance of several variations of the two basic 
filters developed in Chapter III. Descriptions and results 
of the Monte-Carlo simulations are presented in Chapter IV. 

In Chapter V the results of the Monte-Carlo tests reported 
in Chapter IV are used to develop two methods of making the 
filter adaptive to changing situations. Monte-Carlo simula- 
tion results for the adaptive filters are reported in Chap- 
ter VI. 

The final chapter summarizes the results of the various 


Simulations and presents conclusions. 


B. BACKGROUND MATERIAL 

The model used for the system, composed of a target and 
a radar which makes noisy measurements of target range, 
bearing and elevation iS approximated by the linear differ- 


ence equations 


x GRE) x(k) = Tw(k) Cr) 


Z(k) = Hx(k) + v(k) ШІ! 2) 


where 


x(k) is the n-dimensional state vector at time 
š t-kT 


W(k) is the p-dimensional random forcing vector at 
b time t-kT 


Ds the nxn istate transition matrix 
T(k) is the nxp matrix which related the random 


forcing terms to the state vector 


EZ 





Z(k) is the m-dimensional vector of measurements 
taken at time t=kT 


H is the mxn observation matrix 


v(k) is the m-dimensional vector of measurement 
d noise 


T is the sampling period 


is a non-negative integer which defines the 
current sample. 


The assumed random process statistics are summarized by 


the following equations. 


E{v(k)} = E{w(k)} = О tor all k CES) 
E(Y (Y (3)) 7 RG08,, (1.4) 
PQOEUQOw G)) (9 = QC) 85, (1.5) 
Eiv(k)w(j)} = 0 fortfall j,k. (1.6) 


ж? ero 
“ç x. mæ мн” 


E is the expectation operator 

E is the Kronecker delta function 

R(k) is the mxm covariance of measurement noise 
Е matrix 


Q(k) is the nxn state excitation covariance matrix. 


The discrete Kalman filter equations are [5] 


G(k) = P(k/k-1)H [HP(k/k-1)H + R(k)17? GT) 
P(k/k) 7 [I - G(k)H]P(k/k-1) (1.8) 
P(k/k-1) = 9P(k-1/k-1)9  * Q(k-1) (9) 
X(k/k-1) - 9Х(к-1/к-1) (LO) 
S(k/k) 7 R(k/k-1) + G(K)[Z(k) - HX(k/k-1)] ge 


Where the notation (k/k-1) denotes the calculation of a 


quantity at time t-kT based on information up to and 
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including time t-(k-1)T. The previously undefined matrices 


which appear in the Kalman filter equations are defined as 


follows: 


G(k) 
I 


s 


P(k/k) 
P(k/k-1) 


X(k/k) 


X(k/k-1) 


the nxm Kalman filter gain matrix, 
the nxn identity matrix, 


the nxn covariance of estimation error 
matrix, 


the one-sample prediction error covari- 
ance matrix, 


the nx1 state estimate vector, 


the nx1 state prediction vector. 


1.4 











II. THE ROTATED CARTESIAN COORDINATE SYSTEM 


A. MOTIVATION FOR ROTATING THE COORDINATE SYSTEM 
1. Correlation of Measurement Noise 
When radar measurements taken in spherical coordi- 
nates are converted to cartesian coordinates for filter 
processing the measurement noise becomes correlated in the 
process. 


The measurements in spherical coordinates are given 


by 
Z (k) = r(k) + n (K) (2.1) 
Z, (k)  B(k) * n,(k) (ОО) 
Z,(k) = E(k) + n,(k) (2.3) 


where r(k), B(k) and E(k) are true target range, bearing and 
elevation respectively and n,(k), no(K) and n3(k) are the 
corresponding measurement noise components. It is assumed 


that пу (К), п, (к) апа роз (к) are mutually uncorrelated and 


that 
Ein (k)) 7 E(n,(k)) = E(n (&)) = O 
Efn; (k)) — of, E{n?(k)} = of, E(nt(k)) 7 o4 , 


kEe>2077,2853,,. 22. 


[f the coordinate converted measurements are written as 


i 


Zi(k) = x(k) + vi(k) 


Z, (k) 


II 


y(k) + və(Kk) 


Za(k) 


ACK) * va(k) 
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it can be shown that [2] 


E{v(k)} = 0 


Е{у(К)у (К)} = (К) 


where the entries of the covariance matrix of measurement 


error R(k) are 


> 


Rii(k) = r'(k)opcos^B(k)sin?E(k) 


114: 


0:9 
R, (K) 


R,, (k) 


R,, (K) 
R, (K) 
R,, (k) 


R33(k) 


+ r*(k)ogsin*B(k)cos*E(k) 
+ r? (k)ogofsin*B(k)sin*E(k) 
Í 0 cos*B(k) 

2 : 1.2 И 2 
r“ (k)cosB(k)sinB(k)sin E(K) [05 овор] 
+ ѕіпВ(к)соѕВ(к)соѕ?Е(К) [05 -г? (К)ср] 
cosB(k)sinE(k)cosE(k) [o2 -r? (k)o£.] 
К, (К) 
r^(k)giísin'B(k)sin^E(k) 

+ r* (k)o,cos*B(k)cos? E(k) 
+ r* (k)ogocos*B(k)sin*E(k) 
+ 0 sin*B(k)cos* E(k) 


sinB(k)sinE(k)cosE(k) [g2-r? (k)or.] 


R Ck) 


R, , Ck) 


r*(k)0,cos*E(k) + P in Oh) 
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( 2 


(2. 


(2. 


OM 


(2. 


(2. 


22 


(2. 


(2.4) 


СО) 


б) 


7) 


8) 


9) 


tO) 


100) 


12) 


3) 





The GFCS MK 86 filter assumes zero mean uncorrelated noise 
[1], however, equations (2.5) through (2.13) indicate that 
in general the measurement noise is correlated. 

The effect of the correlation is diminished in the 
rotated cartesian coordinate system, since for an inbound 
target the bearing will be zero or near zero. If B(k) = O 


is assumed, 


r*(k)or. sin? E(k) sinE(k)cosE(k)[or 


2 ep2 2 
HO r°(k)on] 


r? (k)o2[ cos’ E(k) 


R(k) = O О 
+topsin*E(k)] 
sinE(k)cosE(k)[o., r? (k)o pcos E(k) 
О 
-r? (k)o2] +o sin E(k) | 
ü J 
(2.14) 


Although the measurement noise is still correlated, the de- 
pendence on target bearing has been eliminated. For the 

case of the low elevation target, where the measurement noise 
is expected to be relatively high due to poor tracking con- 
ditions in the radar, the approximations sinE(k)=0 and 


COSE(k)*1 are valid and 


Gre О О 
r 
R E. 2 2. 
925) 21% r'(k)9g 0 
0 0 ré qo i 15) 


which agrees in form with the uncorrelated noise model. 
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For the fixed cooreiinate system the low elevation assumption 


yields 
r^(k)o5sin^B(k) sinB(k)cosB(k) [o7 
0 
2 2 : u 2 
+0 Cos B(E) r(k)og] 
r^(k)o$cos^B(k) 
R(K) = Riko, 0 
~ Ze? 
+0, sin B(k) 
0 О ro(k)0, 


(2.16) 


which agrees with @he uncorrelated noise model only when 
B(k) = O, 90, 180, or 270 degrees. Equations (2.15) and 
(2.16) indicate that the uncorrelated measurement noise as- 
sumption is likely to be considerably more accurate in the 
rotated nern than it is in the fixed cartesian 
coordinate system. 
2. Modeling of Random Forcing Inputs 

The relative magnitudes of the vector components of 
forcing inputs in fixed cartesian coordinates depend on the 
target's orientati@n within the coordinate system. 

The track generating program described in Appendix 
C has been used ta demonstrate the change in the relative 
magnitude of the x and y acceleration components of a target 
maneuvering in the xy plane as the bearing of the track's 
axis of symmetry is changed as illustrated in Fig. 2 for two 
bearings. Inbound trajectories with seven-g maneuvers in 
the xy plane were generated using the track generating pro- 
gram in ten separate runs with each run differing only in 
the bearing angle wf the axis of symmetry of the trajectory. 
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Table I shows the relationship between the maximum 


, and the maximum 


acceleration in the x-direction, (| m 


acceleration in the y-direction, |y | for each run. 


max' 


TABLE I 


COMPARISON OF MAXIMUM ACCELERATIONS 


BEARING OF AXIS 


= OF SYMMETRY (deg) IX | ax (E? APER, 
1 О 0.6 7.0 
2 10 1.4 6.9 
3 20 2.6 6.7 
ч. 30 3.6 ОТЕ 
o 40 4.7 9.4 
6 o0 5 4.6 
F 60 6.2 э 
& 70 6.8 2.0 
3 80 6.9 1.4 
10 90 7.0 0.7 


The large variation in the relative magnitudes of 
the acceleration components indicates that much more accu- 
rate apriori estimates of random forcing inputs are possible 
if all targets approach the origin from the same bearing. 
Therefore a considerable advantage in filter design might be 
gained by rotating the coordinate system so that the x-axis 
is aligned with the initial detection bearing of each new 


target. 
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B. COORDINATE SYSTEM ROTATION 

The radar measurements of target position are range, 
r(k), bearing, B(k), and elevation, E(k) (to simplify the 
notation in the subsequent discussion r(k), B(k), and E(k) 
are assumed to include the measurement noise). To convert 


the measurements to cartesian coordinates the equations are 


x(k) = r(k)cosB(k)cosE(k) (2.17) 
y(k) » r(k)sinB(k)cosE(k) (2.18) 
Z(k) = r(k)sinE(k) (2.19) 


where once again for simplicity of notation x(k), y(k) and 
Z(k) are assumed to include measurement noise. 
If the x-axis is aligned with the initial bearing of the 


target, equations (2.17) through (2.19) become 


X'(k) = r(k)cosAB(k)cosE(k) (2.20) 

y'(k) = r(k)sinAB(k)cosE(k) (2.21) 

z'(k) = r(k)sinE(k) (2.22) 
where 

AB(k) = B(k) - B(0) | (2.23) 


as depicted in Fig. 3. 

In order to maintain the orientation of the rotated co- 
ordinate system it is necessary to store the initial bearing, 
B(O), for use in the coordinate conversion process as indi- 


cated in equation (2.28). 
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FIGURE 3  ROTATED COORDINATE SYSTEM GEOMETRY IN RELATION TO 
FIXED CARTESIAN COORDINATES 
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III. FILTER DEVELOPMENT 


A. THE CONSTANT-VELOCITY (1/S?) MODEL 

The first filter considered is based on a model which 
assumes that the target maintains constant velocity in all 
three coordinate directions. 

The state transition matrix, Ф, апа г, the matrix which 
relates the random forcing inputs to the state vector are 


[1] 


TET 
O 1 he 
EU I 
o(T) = O | | О (3.1) 
LU 9 1; =e 
"MU 1gT 
U ) ` i 
Lei ejoa 
and 
2 | 
T'/2 | О i O0 
T "V NT, 
>= - - I- —— |— — — 
| 2 ! 
^ ~ | T | 2 
| 
шу m — 
EE 
р | 


for the case oí the constant-velocity system model with 
sampling time T. 

The covariance of measurement noise matrix R(k), for the 
system is given in Chapter II, equation (2.14). If the fil- 
ter gains are to be calculated off-line, a model must be 


chosen for R(k). This necessitates that values be chosen 
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for target range, target elevation, and the standard devia- 
tions of measurement noise. Range and elevation values were 


arbitrarily chosen as 


r(k) 


10 kyd (3.3) 
E(k) 


O deg (3.4) 
the measurement noise standard deviations used were 


0. = a yd (3.59 


E = E - 0.3 mrad (3.6) 


These values of On? Op: and соь are considered to be typical 
of good tracking radar systems.  Substituting the values 
given in equations (3.3) through (3.6) into equation (2.14) 
yields 
"= c 
9 O О 
кокус. O 9 О (507) 


О О 9 


The system model described above is highly simplified. 
The model allows only for targets which fly at constant ve- 
locity. It is not expected that any target which the system 
might be used against would behave so well. Almost certain- 
ly there will be pilot initiated accelerations, wind buffet- 
ing, etc. which are not accounted for in the model. It is 
the function of the state excitation covariance matrix Q tO 
compensate for the inadequacies of the simplified model. 

The starting point for developing the state excitation 
covariance matrix is the random cie input vector, W, 


represented by 





w= a (3.8) 


where a a and a, are accelerations in the x, y, and 
z-coordinate directions respectively. Since any coupling 
between the accelerations depends on the type of random 


forcing experienced, it is assumed for this model that the 


random forcing is decoupled, i.e. 
Eta, a) = Е(а,а,) = A же = 0 
The state excitation matrix Q is defined in Chapter I as 


(3-2) 


If the matrix multiplication is performed with ! and w as 
previously defined, and the expectation operation períormed 


on the product, the result is 


E 
X X 























| | 
| 
4 2 : | 
О О 
w manu ~ | ^: 
= M | 
2 | | 
— — —=—————— —-———— 
| а^т° ae | 
K 
| 
| aT? айта | 
| 2 | 
айке С. Фа. --- - - - | -- - - - - 
| | a 2 p" 2913 
| | 7, 
О | О | 
| | 
I | 


a an P 
7, 2, 
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The values assumed for a DE and a, for off-line cai- 
culation of the filter gains have a significant effect on 
the Hitter s performance. Much OTe thnessimu lation work re- 
ported in Chapter IV is devoted to the study of the effects 
of different assumed values of a. == and а, on the filter's 
performance. 

The Kalman filter requires that an apriori estimate of 
the system state vector be made. The initial state estimate 
vecrcorsucdEnaartlossrmulation runseis?7x(0/—199— O, based on 
the reasoning that in a fixed cartesian coordinate system 
this is the mean of all possible initial state vectors. 
Since the system being simulated assumes a rotated cartesian 
coordinate system as described in Chapter II, it might be 
argued that спе x coordinatc position could be estimated as 
some anticipated detection range, and the x-coordinate ve- 
locity as a typical aircraft speed in the negative x-direc- 
tion. However, since the filter converges rapidly as it is 
presently initialized, it is felt that to revise the initial- 
ization procedure would only serve to complicate the model 
and would provide no foreseeable benefit. 

What is believed to be the simplest approach is taken 
again in choosing the initial covariance of estimation error 
matrix. It is desired that the initial position measure- 
ments be used as the initial position state estimates in 
order to Nave the filter track as quickly as possible. This 
end is accomplished quite simply by making the diagonal ele- 


memts Oi the initial covarimnce of estimaition error matrix 





very large (e.g. 10?). The effect of doing this can be 


seen by examining the Kalman filter gain equation 
» 1 T =m] 
G(k) = P(k/k-1)H [HP(k/k-1)H +R(k)] CS du 
and considering the scalar case 


G(k) P(k/k-1)H[HP(k/k-1)H+R(k)]7? 


. | P(k/k-1)H (3.12) 
H?P(k/k-1)+R(k) 


where since H is the observation matrix the scalar equiva- 
lent would simply be unity. If P(k/k-1) is made very large 


with respect to R(k) for k-0, then we have 


P(0/-1) 


EXOD) Sk Ser 1) In (0) — 


1 (СО; 


which makes the first estimate equal the first measurement -- 
a logical result in the presence of a high degree of uncer- 
tainty. 

The initial covariance of estimation error matrix used 


íor all simulations is 


10? 
10? 


TO 


P(0/-1) = E 


пох (Sov) 
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B. THE CORRELATED RANDOM ACCELERATION MODEL 

One approach used to simulate the effects of target 
motion involves introducing correlated random acceleration 
buts to a 1/S? model. The correlated random accelerations 
are generated by an input of white noise into a coloring 
filter having a transfer function of a/(Sta). The value of 
a can be selected to vary the time constant of the maneuver. 
This results in an overall transfer function for the model 


of target motion of 


T(S) = (3. B9) 


EN __ 
s*(S+a) 


The 9x1 system state vector assumed for this model of 
target motion is 


Гесу) 

x(t) 

alt) 

y(t) 

x = Ту» 
y(t) 

z(t)! 

z(t) 

Z(t) | (3.16) 





McKinley [4] describes this model in detail and gives 
the state difference equations (for one dimension) from 


which the state transition matrix 
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1 m i = + ШЕ, 
ФС) = [о а 21-е 
n о “Тут (3.17) 
and 
г,(Т) = | т - i ES] 
шыс е! (3.18) 


can be extracted. 
The matrices 2.(Т) апа r CT) apply to a three-state sys- 
tem in one coordinate direction. For the full, nine-state 


three-coordinate system, the 9x9 matrix %@(T) is 


OO 0 0 
MOET Soe aT) 7B 
О О OG) Eco) 


where $ (T), d and ó (T) are each equal to the 3x3 ma- 


Ir lx ó, (T), and the 9x3 matrix [(T) is 


n 0 0 
P го) о 
0 0 Г СТ) (3.20) 


Le. 
where r. CT), LS and r. CT) are each equal to the 3x1 vec- 
тот P (T). 
Jn order to make the derivation of the Q matrix less 


ООН Л СП tO follow the following notation ist adopted. 
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үй 0 0 
5 О О 
de О О 
О jr О 
en) 0 ve 0 
0 Үз 0 
О О \ 
О О Y. 
0 0 Y. ӘЛІП 
where р 
it T 1 -aT 
y; el deae (3.22) 
PERLE -= ІЗЕТ ы , and (3.23) 
Yo Qd Ges) 


The forcing input vector w and the state excitation co- 
variance matrix Q are defined as in equations (3.8) and (3.9), 
respectively. If the matrix multiplication is done and the 


expectation operation performed on the product, the result 


is 
ах 4 0 0 
жр 22d ° 
0 0 а> а au 25) 
where 
Y [Үл KGY, 
О Y. M 
ЛӘ Yates y (3.26) 
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The assumed covariance of measurement noise matrix, R, 
is the same as given in equation (3.7) for the constant- 
velocity system. 

The constant a is inversely proportional to the assumed 
maneuver time constant. For example, the value of a=0.5 
corresponds to a time constant of two seconds which is con- 
sidered to be typical for an aircraft evasive maneuver. 
Different time constants could conceivably be applied to 
each coordinate direction, however, in this thesis one value 
of a is applied throughout the system model. 

The a/S*(Sta) model was simulated using values of a=0.1, 
0.3, 0.5 and 1.0. Based on the results of the simulations -- 
reported in-Chapter IV -- a=0.5 was chosen as the value pro- 
viding tne best overall performance in the limited numher of 
tests conducted. This does not suggest that a=0.5 is an 


optimum choice. 





IV. MONTE-CARLO SIMULATION 


` 


The Monte-Carlo simulation program described in Appendix 
A was used to evaluate variations of the two basic models 


described in Chapter III. 


A. TEST TRAJECTORIES 

A set of four trajectories was chosen as a standard for 
use in simulating all models in order to facilitate compari- 
son of the simulation results. All of the trajectories are 
products of the track generating program described in Appen- 
dix C. Each track is comprised of 60 state trajectory data 
samples with a sampling rate of four hertz. 

The tracks are developed so that the initial target 
bearing coincides with the x axis to simulate the rotated 
coordinate system. The addition of measurement noise will 
cause the initial bearing to not coincide exactly with the 
X axis, however, no significant error is expected. 

1. The Short-Range, Constant-Velocity Trajectory 

The short-range, constant-velocity trajectory is a 
Stragehesseinbound trdjectory, beginning at a range of 3,3350 
yards at 30 degrees elevation and proceeding toward the ori- 
РЕПИ а frange ratefof 230 yd/sftofa Tınalgrange of 460 
yards. 

Pee Lites Lone—Rangse weconstant—-Velocity Trajectory 

The long-range, constant-velocity trajectory is a 


straight inbound track beginning at a range of 9,225 yards 





at 30 degrees elevation and proceeding toward the origin 
with a range rate of 230 yd/s to a final range of 5,817 
yards. 
3. The Short-Range Trajectory With Maneuver 
The short-range trajectory with a maneuver, like the 
short-range, constant-velocity track, begins at a range of 
3,850 yards at 30 degrees elevation and proceeds toward the 
origin with a target speed of 230 yd/s. Commencing at a 
range of 3,500 yards the target undergoes a two-g, sin?’ 
maneuver of 13.5 seconds duration. The maneuver is completed 
at the final sample range of 470 yards. 
4. The Long-Range Trajectory With Maneuver 
The long-range trajectory with a maneuver begins at 


"e m m y „= v 


= СУИ У = ш Mm 
a Launga Ul „шшс y L 


aras av 30 Gegrees elevation. The target 
proceeds toward the origin with a speed of 230 yd/s. Com- 
mencing at a range of 8,900 yards the target undergoes a 
two-g, sin* maneuver of 13.5 seconds duration. The maneuver 
is completed at the final sample range of 5,853 yards. 
5. Other Trajectories 

Representative variations of the two basic models 
were simulated using various trajectories in addition to the 
four standard trajectories described above in order to in- 
Sure that the Characteristics shown in thel simulations were 
valid for trajectories other than the standard set. 

The additional trajectories used include a constant- 


velocity erossing trajectory which has.its initial point on 


the x-axis (simulating the rotated coordinate system) at a 
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range of 3,500 yards and proceeds to diverge from the x-axis 
at an angle of ten degrees as range decreases, and inbound 
trajectories with two-g and 2.5-g sin? maneuvers of ten sec- 


onds duration occurring at various ranges. 


B. ACCELERATIONS ASSUMED 

The variations of the two basic system models differ 
primarily in the values of accelerations assumed for the 
components of the forcing input vector, w(k), which is used 
in calculating the state excitation covariance matrix QCk). 
As indicated in equation (3.10) for the constant-velocity 
model, and equation (3.25) for the correlated random accel- 
eration model, larger assumed values of acceleration result 
in larger corresponding entries in the Q matrix. Solution 
of equations (1.7) through (1.9) shows that a consequence of 
larger Q matrix entries is larger entries in the gain matrix 
G(k). This results in the filter being more responsive to 
target maneuvers, but it is also more susceptible to estima- 
tion error caused by measurement noise. The assumption of 
different values for the components of w(k) is done in an 
effort to determine the Q matrix which best serves to com- 
pensate for modeling inadequacies as discussed in Chapter 
III. 

The values of acceleration assumed are chosen as being 
representative of manned aircraft accelerations in the ro- 
tated cartesian coordinate system. It is predicted that an 
attacking aircraft during evasive maneuvering would have ac- 


Scrertcions іп they and z-directions Or three g*s.or less 
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and accelerations in the x-direction of less than one g. 
The values used for the simulation runs are listed in Tables 
II through IV along with the corresponding simulation pre- 


diction accuracy statistics. 


C. SIMULATION RESULTS 

The main criterion chosen for evaluation of model per- 
formance is the position-state prediction accuracy as dis- 
cussed in Appendix B, although the means and covariances of 
estimation error are also calculated. The predictions for 
each model are calculated using the state transition matrix 
for the model as developed in Chapter III. 

The statistics tabulated for comparison are the percent- 


age of predictions per run which were in error by less than 


24 feet -- labeled "hit %'' -- and the percentage of predic- 
tions per run in error by less than 48 feet -- labeled 
| care 9". 


The results of the various constant-velocity model simu- 
lations are listed in Table II. All statistics listed are 
products of 100 member Monte-Carlo ensembles. 

The correlated-random-acceleration model was tested for 
several values of a as well as various assumed accelerations. 
The simulation results for a=0.3 and a=0.5 are listed in 
Tables III and IV respectively. All other values of a 
tested provided inferior performance and the results are not 
tabulated. 

The simulation results listed in Tables II through IV 


are considered to be representativo of all prediction 
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accuracy statistics obtained, with the exception that Table 
II does not present statistics indicating the effect on fil- 
ter performance when the assumed x-direction acceleration 

a, is varied. It was determined early in the simulation 
phase that the constant-velocity model performs best against 
the test trajectories used when a, is approximately O.1 g, 
therefore in the interest of computer time conservation oniy 
this value of a, was used for the final set of simulation 
runs. 

Also in order to conserve computer time the correlated- 
random-acceleration model with a=0.5 was not tested against 
the long-range trajectory with a maneuver. This is consid- 
ered a justifiable omission Since no appreciable improvement 
or further degradation of prediction accuracy in this area 
is anticipated. 

In Table V the five best performing filters against each 
Lest trajJectory are ranked according to the prediction accu- 
racy statistics shown in Table II through IV. The "hit" and 
"Scare'" percentages for each ranked filter are included with 
the model description. For the long-range, maneuvering tar- 
get only two of the estimators tested produced "scare" per- 
centages of one percent -- the filter ranked highest is 
chosen because it produced smaller average miss distances 
than the second and third ranked filters for 17 of 21 pre- 
dictions per run even though it scored only two "scares" in 


the 100-member ensemble. 
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D. COMPARISON OF PREDICTION ACCURACY AND ESTIMATION ERROR 

STATISTICS 

An effort to correlate the prediction accuracy statis- 
tics with the estimation error means and variances for the 
sample estimate from which the prediction is made indicates 
that for constant-velocity targets, where the estimation er- 
ror means are small for all the models, the model with the 
smaller variances of estimation error -- especially for the 
velocity state estimates -- can be expected to produce more 
consistently accurate predictions. For example, a compari- 
son of the statistics for two 1/S* models, one with forcing 
inputs assumed to be zero and the other with assumed values 
of O.1 g in the x-direction and two g's in the y and z-di- 
rections, shows that both filters in general have small 
means of estimation error. The estimate at sample number 18, 
x(18/18), upon which the prediction of target position at 
the time of sample number 53, &(53/18), is based in the long- 
range, constant-velocity trajectory, has small means of es- 
timation error for both filters (less than one yard for 
position states and one yd/s for velocity states). The var- 
iances of y and z-direction velocity estimation error for 
the model which assumes two-g acceleration in the y and z 
directions are much larger (an order of magnitude larger for 
the y and z-direction velocity states). This model scores 
18 "scares" at sample number 53 over the 100 run ensemble 


while the model which assumes zero forcing inputs scores 100 
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"scares". The trend indicated in this example holds for all 
constant-velocity trajectory statistics tested. 

A similar analysis of the statistics generated for short- 
range maneuvering targets, where the means of estimation er- 
ror are in some cases more than 20 yards for the position 
states and more than 30 yd/s for the velocity states, indi- 
cates that the better prediction accuracy statistics are 
produced by the model with the lower means of estimation er- 
ror. The variance of estimation error does not appear to be 
a significant factor at short range in any of the models 
tested. Notice that the mean estimation errors are much 
larger than the measurement errors, this is due to the tar- 
get maneuver occurring after the constant-velocity filter 
has reached steady state. 

At long ranges, any "hits" scored against a maneuvering 
target are considered to be either purely by chance or due 
to the type of maneuver which the target executes (i.e. the 
target may pass close to the predicted position even though 
it has maneuvered radically during the projectile time-of- 
flight). 

Thefaverage prediction miss distancelis Considered to be 
the best prediction accuracy statistice for use in evaluating 
a filter's performance against a long-range, maneuvering 
target. The model ranked first against long-range maneuver- 
ing targets in Table V has larger means and smaller vari- 
ances of estimation error than the second and third ranked 


models. These statistics indicate that the estimator is not 
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tracking the target accurately, but is generating predictions 
of future target position in the right general area. This 

is suspected to be largely because the target maneuvers 
about an axis of symmetry in the test trajectory rather than 


at random. 
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V. ADAPTIVE FILTERS 


A.' THE SWITCH-ON-RANGE ADAPTIVE FILTER 

The simpler of the two adaptive filtering methods pro- 
posed is the switch-on-range adaptive filter. This scheme 
is based on the idea that since the measurement noise in 
cartesian E ES S proportional to the target range, 
relatively high filter gains can be used on a short-range 
tzrget without significant degradation of the estimate qual- 
ity for a constant-velocity target and with significant im- 
ME UE Of the estimate quality for a maneuvering target. 

Upon target detection the estimator uses a gain schedule 
which is based on a small forcing inputs -- such as O.1 g in 
each coordinate direction. This "low Q" gain schedule is 
continued until the target range is less than a predeter- 
mined threshold. Upon reaching the switching threshold 
range the estimator switches to a steady-state gain matrix 
which is based on relatively large forcing inputs ("high Q"). 

Phesl/S° model is used for both the “Jew о“ апа "high б" 
filters. This model was chosen because of its superior per- 
formance against constant-velocity targets, especially at 
long range, as indicated in the Monte-Carlo simulation re- 
suits of Chapter IV. 

The "high Ө" filter uses a steady-state gain matrix be- 
cause it is assumed that a good estimate is available at the 


time the switching occurs. 
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The switching threshold range for simulation purposes is 
chosen as 3,333 yards. The simulation results of Chapter IV 
indicate that the "high Q" filters based on the 1/S? model 
with assumed accelerations of one to three g's in the y and 
z-directions perform relatively well within this range 
against either maneuvering or constant-velocity targets. 

The "high Q" gain matrix used in the simulation is the 
steady-state gain matrix of the 1/S^ model with assumed ac- 
celerations of 0.1 g in the x-direction and three g's in the 
y and z-directions. The gains are computed from the model 
developed in Chapter III. 

A FORTRAN listing of the switch-on-range adaptive filter 


subroutine is included in Appendix D. 


B. SWITCH-ON-RESIDUALS ADAPTIVE FILTER #1 

The Switch-on-residuals adaptive filter is more sophis- 
ticated and more complex than the switch-on-range filter 
ce ly described. 

The first version of the switch-on-residuals adaptive 
filter is designed to track with a "low Q" gain schedule 
based on the constant-velocity model until a component of 


the filter residual vector 
p(k) = Z(k)-Hx(k/k-1) (oa 
exceeds a calculated gating level, 


а (К) = 2/P,. Ck/k)HR, (Ck) , i=1,2,3, ( 


сл 
ho 
м2 


on two successive samples. 
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P.,(k/k) is the theoretical variance of position esti- 
mate error and R. ; (K) js the theoretical variance of measure- 
memt noise in one coordinate direction. 

The values used for the theoretical variance of position 
estimate error are taken from the steady-state theoretical 
covariance of estimation error matrix for the model from 
which the “high Q" gains are calculated. 

The values used for the theoretical variance of measure- 
ment noise in the gating level calculation are calculated 
on-line from the diagonal terms of the covariance of meas- 
urement noise matrix given in equation (2.14) and vary with 
the target location. 

If the occurrence of a residual value greater than the 
gating level at time i is defined as event А; and the prob- 
abiiity of this occurring on the next sample is event Asya 
A; anq А. are independent events. The probability of two 
successive samples with residuals greater than the gating 


lewel is 
Р(А А, 41? = PA IPLA 44) à 65.39 


If only the R. (K) term of the equation (5.2) is consid- 
ered, then the gating level is equal to twice the standard 
deviation of measurement noise. For a normal distribution 
the probability of an event differing from the mean by more 


than two standard deviations is givon as 


P(A.) = PY = 0.0228 


A544) 


46 





Substitution of this value into equation (5.3) yields a 


probability of 


4 P(A,A ) = (0.0228)? = 5.20 x 10" 


dk 
which is the probability that the residual is greater than 
the gating level due to estimation error or measurement 
noise when no target maneuver is experienced. 

The filter is лери лет; adaptive in each coordinate 
direction. For example, when p,(k) is greater than £ (K) 
for two successive samples, the x-direction states are re- 
initialized by accepting the latest measurement as the posi- 
tion estimate and beginning at k=0 with a "high Q" gain 
schedule based on the a/S*(Sta) system model with values of 

and assumed forcing inputs chosen to provide good perform- 
ance against a EN ins target. 

The y-direction and z-direction state estimates continue 
to be based on the constant-velocity (1/S?) model. The 
a/S?(S*g) filter gains continue to be used in making the 
x-direction estimates until three successive samples have 
residuals less than one-half the gating level. 


Then the inequality 
p(k) < VP,,CkE/k)*R,  (K) (5.4) 


Beats Tea Lor threefsúuccessive Samples the filter switches 
back to: the "low Q" gain schedule assuming the same time in- 
dex as would have been used for the current sample if no 


adaptive switching had occurred. 
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For the Monte-Carlo simulation of the switch-on-residu- 
als filter, the "low Q" gains are based on the 1/S? model 
with assumed forcing inputs of 0.1 g in each coordinate di- 
rection since it is proposed tbat those gains will be used 
exclusively against constant-velocity targets. This choice 
is based on the simulation results listed in Chapter IV. 

Also based on the simulation results listed in Chapter 
IV, the model chosen for calculation of the "high Q" gains 
is the a/S*(Sta) model with a-0.5 and assumed forcing inputs 
of 0.1 g in the x-direction and one g in the y and z-direc- 


tiens. 


C.  SWITCH-ON-RESIDUALS ADAPTIVE FILTER #2 

The second version of the switcb-on-residuals adaptive 
fitter differs irom the first in the manner in which the 
Switching gates are calculated. 


In this version the gating level is calculated as 
* = " o . < ] = . 
2. (K) SYP. (К/К)+Е.. (К) A (5.5) 


an the switch is made to the "high Q" gain schedule when 
the gating level iS exceeded by one sample residual. 

The switch back to the "low Q" gain schedule occurs when 
the residual is less than 0.88 & (K) for two successive sam- 
ples (the coefficient 0.88 is an arbitrary choice for a val- 
ue greater than 0.5 but less than one). In this version 
when the filter switches back to the "low Q" gains the index 
fion the "high Q" gain schedule is not reset immediately but 


1s delayed one sample so that if the residual exceeds the 
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gating level on the first sample after the filter returns to 
the "low Q" gain schedule, the filter does not reinitialize 
by accepting the position measurement as the estimate but 
simply switches back to the "high Ө" gain schedule using the 
same time index as would have been used if the switch to the 
"low Q' gain schedule had not taken place. 

The gain schedules used in this filter are the same as 
those used in the first version of the switch-on-residuals 


adaptive filter described in the previous section. 


D. THE COMPROMISE FILTER 

The compromise filter is an attempt to pick one gain 
schedule which provides adequate performance for all types 
of trajectories without adaptive switching. The filter 
simulation results listed in Chapter IV indicate no clear 
пош е ог thes compromise filter. 

Because of its good performance against constant-veloc- 
ity targets at all ranges and its acceptable performance 
against short-range maneuvering targets, the 1/S* model 
which assumes accelerations of 0.1 с in the x-direction and 
De O2 Sm леу апа z-directions is ehosen for calculation of 


the compromise filter gains. 


E. SUMMARY OF ADAPTIVE FILTER DEVELOPMENT 
Four adaptive filters have been developed for evaluation 
by Monte-Carlo simulation, Table VI provides a summary of 


their key features. 
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VI. ADAPTIVE FILTER SIMULATION 


"The Monte-Carlo simulation program described in Appendix 
" A is used in conjunction with the adaptive filter subrou- 
tines listed in Appendix D to evaluate the performance of 
the three adaptive filter schemes and the compromise non- 


adaptive filter proposed in Chapter V. 


A. TEST TRAJECTORIES 
A set of six test trajectories was used for the simula- 
tion of each of the adaptive filter configurations and the 
compromise filter. Each test trajectory is a product of the 
track generating program described in Appendix C and is com- 
prised of 160 state trajectory data samples with a four- 
hertz sampling rate. All six of the test trajectories 
represent a target proceeding inbound from a range of 6,000 
yards and an elevation of 20 degrees at a speed of 200 yd/s 
with a dive angle of 20 degrees. The trajectories differ in 
the target maneuvers incorporated: 
Trajectory No. 1 - No maneuver 
Trajectory No. 2 - À two-g maneuver of 13.5 seconds 
duration commencing at a range of 
4,000 yards. 
Trajectory No. 3 - A two-g maneuver of 13.5 seconds 
duration commencing at a range of 
2 TADO ER 


Trajectory No. 4 


A two-g maneuver of 13.5 seconds 
duration commencing at a range of 
3 333 yards. 
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Trajectory No. 5 - A 2.5-g maneuver of 10 seconds dura- 
tion commencing at a range of 2,700 
yards. 

Trajectory No. 6 - A 2.5-g maneuver of 10 seconds dura- 
tion commencing at a range of 4,000 
yards. 

B. SIMULATION RESULTS 

The prediction accuracy statistics produced by the Monte- 
Carlo simulation of the switch-on-range and switch-on-resid- 
uals adaptive filters and the compromise constant-velocity 
filter are listed in Table VII. The statistics presented 
are products of 80 member ensembles. 

The adaptive filter predictions are generated using the 
constant-velocity model state transition matrix for all con- 
figurations. A small number of simulation runs made using 


the correlated-random-acceleration model state transition 


matrix produced slightly lower "hit" and "scare" percentages. 
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VII. SUMMARY AND CONCLUSIONS 


A. TEST RESULT SUMMARY 
1. Simulation Result Summary for the Two Basic Models 

The simulation results presented in Chapter IV show 
that the constant-velocity (1/8?) model with relatively 
small values of assumed forcing inputs provides the better 
performance against constant-velocity targets at both short 
range and long range. By increasing the value assumed for 
the forcing inputs to the 1/S* model its performance against 
maneuvering targets can be improved but there is a corre- 
sponding loss in the performance for constant-velocity tar- 
gets, particularly at long range where the transverse 
measurement noise components are reiatively large. 

The correiated-random-acceleration model provides 
better performance against maneuvering targets and reason- 
ably good performance against constant-velocity targets at 
short range, but its performance against long-range, con- 
stant-velocity targets, as measured by the prediction accu- 
racy statistics listed im Chapter IV Ts imterlor Lo the 
"low Q" 1/S? models. 

None of the models tested provided good prediction 
accuracy against the long-range maneuvering target. This is 
to be expected because of the long time of flight over which 
the prediction must be made. Even a perfect estimator would 
surely provide inaccurate predictions when the predictions 


must be made ten seconds ahead on a target which maneuvers 
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in the interim. Average miss distance was chosen as the 
most appropriate statistic for evaluating a filter's per- 
formance against a long-range, maneuvering target. 

2. Adaptive Filter Simulation Summary 

The simulation results for the adaptive filters and 
the compromise filter listed in Chapter VI indicate that the 
switch-on-range adaptive filter provides the best overall 
performance against the six targets trajectories included in 
the test. Although the performance of the second version of 
the switch-on-residuals adaptive filter is comparable. 

It was anticipated that the more sophisticated 
switch-on-residuals adaptive filter would provide the better 
performance, however, even the second version of the switch- 
on-residuals filter, which was partially developed by "cut 
and try" modifications to the original version, does not 
better the performance of the simpler, switch-on-range 
adaptive filter. The problem encountered with switch-on- 
residual filter seems to be that measurement noise causes 
the filter to switch to the "high Q' gain schedule when 
tracking long range constant-velocity targets. 

The compromise filter, the simplest approach, pro- 
vides fair performance -- better than the original version 
of the switch-on-residuals adaptive filter -- and serves as 
a standard for comparison with the adaptive filter perform- 


once. 
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B. SUGGESTIONS FOR FURTHER INVESTIGATION 

One of the most obvious extensions of the simulation 
work is to incorporate the rotated coordinate system algo- 
rithm developed in Chapter II into the Monte-Carlo program 
for complete testing. 

The switch-on-residuals adaptive filter should receive 
more attention. This is believed to be potentially the bet- 
ter of the adaptive filtering schemes but the most effective 
combination of gating levels and switching rules has not 
been found. 

The switch-on-range adaptive filter shouid be simulated 
using other values for the steady state "high Q'" gain ma- 
trix. A prime candidate is the steady-state gain matrix of 
the correlated-random-acceleration model with a=0.5 and 
assumed forcing inputs of O.1 g in the x-direction and one g 
in the y and z directions. 

The various models should also be simulated using other 
types of trajectories including a proportional navigation 
trajectory representing an anti-ship missile attack. 

More effort should be given to developing a definite 
quantitative correlation between prediction accuracy and 
estimation error statistics as discussed briefly in Chapter 
IV. A thorough investigation of relationships between the 
two sets of statistics may provide useful information for 


use in filter design. 





C. CONCLUSIONS 

The rotated cartesian coordinate system has been shown 
to have potential advantages over the TR EA Las ian coor- 
dinate system both in diminished correlation of measurement 
noise and in selection of the components of assumed forcing 
inputs for the system model. 

The prediction accuracy statistics are concise and valid 
measures of the system's performance in most cases, however 
the mean and covariance of estimation error are valuable 
statistics to have available to support the prediction accu- 
racy statistics and to serve as measures of filter perform- 
ance in cases such as the long range maneuvering target, 
where the prediction accuracy measurement provides little 
information. 

The constant-velocity model with small values of assumed 
forcing inputs has been shown to provide the best perform- 
ance against constant-velocity aircraft while the correlated- 
random-acceleration model with the proper choice of a to 
match the maneuver time constant provides the best perform- 
ance against maneuvering targets. 

Finally, the relatively simple switch-on-range adaptive 
filter and the more complex switch-on-residuals adaptive 
filter have been shown to provide improved performance over 


the compromise non-adaptive filter. 
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APPENDIX A 


THE MONTE-CARLO SIMULATION PROGRAM 


A.1 Program Description 


The general purpose Monte-Carlo simulation program de- 
veloped by Prof. D. E. Kirk of the Naval Postgraduate School 
is quite flexible in its basic form and offers several 
options to the user. The various options and many of the 
input data requirements are explained by comments early in 


the program listing. 





A.2 The Monte-Carlo Simulation Program Listin 
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A.3 


Flowchart of the Monte-Carlo Simulation Program 


N,M, IN, NSAM 
NENS 








IC,IFLR,IEST 
TURK ee nA te. 
IQ,ITRO, IMISS 










IGPLT, IPHPLT 
ІМІРІТ,ІЗШРІЛ 
ISVPLT 
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Iris IH In 
IPKKM1, IGAM, IC 
ISIGV, ISIGW, [XZ 







CALL OVFLOW 
and set kernals 
for SNORM 

IW, IV, IXZ 












Describe the 
run according 
to the options 
selected, 






Indicate the 
input data 
called for 


Indicate the 


output data 
called for. 





'N,M, IN,NSAM, 
NENS. 
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PHI(N,N) 
CALL MREAD 


PHIS(I,J)- 
РН%(1,7) 


PHIX IJI) 
CALL MWRITE 


al? 


H(M,N) 
CALL MREAD 


H(I,J) 








yes 





R(M,N) 
CALL MREAD 








R(M,N) 
CALL MWRITE 








yes 


lI 
= 


IQ 


COVW(IN,IN) 
CALL MREAD 





yes - COVW(IN, IN) 
CALL MWRITE 





Q(N,N) 
CALL MREAD 
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yes 


GAMMA.(N , M) 
CALL MREAD 


no 


GAMMAS (1, J) 
= GAMMA(I,J) 


GAMMA (N,M) 
CALL MWRITE 





[| ^ 


yes 


PKKM1 (N,N) 


ES CALL MREAD 


PKKM1 (N,N) 
CALL MWRITE 





ДА 
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yes 


SIGV(M) 


по E CALL VREAD 





es 
ISIGW= 0 2 


/ V 
` roc | 
SIGW(IN) 


CALL VREAD 


no 


SIGW(IN) 
CALL VWRITE 





XS GE, 1) 
ae estate 


sample 








mo 


ү 


yes 


=i 


no 





XHATZ (N) 
CALL VWRITE 







EN 
XZMEAN (N) 
CALL VREAD 


XZMEAN (N) 
CALL VWRITE 










SIGXZ (N) 









yes 











Read in a 
complete track 
from data cards 
and store in XS 


XS (N, NSAM) 


m iha fin | 


` — _ - 


с 
and last track 
sample points. 


— | 










Form the 
Identity matrix 


EI(N,N) 
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Generate a 
complete track 
and store in XS 





CALL TRACK 





CALL QMAT 


Calculate Q 


using Gammas and 
COVW 


Q(N,N) 
CALL MWRITE 





es no 
2 ТСЕ 


no 






GKS (N,M, NSAM) 
Read filter 
gain schedule 










CALL GAIN 
Generate and 
store PKKS and 
GKS 








Begin 
Monte-Carlo 


loop 





yes 


CALL XZERO 
Generate the 
initial state 
vector 









CALL MEAS 
Form noisy 
measurements 


from state 
values. 
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es 







CALL QON 
Compute Q 


T 


on-line 


yes 


IFLR- 0 


CALL RON 


Compute R 


on-line 


CALL GAIN 
Compute Filter 


gains on-line 















CACL"ESTIM 
Calculate x(k/«) 
and X(k+1)/k) 


CALL STAT 


Compute estima- 


tion error 
Statistics 





K= NSAM TSS 
| no 
es 
no 













CALL TRACK 


Compute Xx (k+1) 


END 
MONTE-CARLO 
LOOP 


O ITER < ENS 
— 


WITER=ENS 
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| Compute the mean 


trajectory 


XM(N,NSAM) 





Compute mean 
estimation error 


ERR(N,NSAM) 


Compute variance 
of estimation 
error 


VAR (N,NSAM) 


no 











Compute off-diag- 
onal covariance of 
estimation error 


IZJ 
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es 


no 


no 





ESSE, ERES 
VAR(I,I,K) 
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= STOP 
yes 
yes 
CALL PLOTP 
k. GKS vs K 
ІТНУРІ=1 \-78 
- CALL PLOTP 
PKKS vs K 
we T 
yes 
no 
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CALL PLOTP 
ERR vs K 





Ж 


ISVPLI- 1 


no 


yes 


CALL PLOTP 
VAR vs K 





STOP 
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A.4 Alphabetical Listing of Key Variables Used in the 


Monte-Carlo Program 


AMISS(NSAM) 


COVW(M,M) 


EI(N,N) 
ERR(N,NSAM) 


G(N,M) 


GAMMA (N,N) 


GAMMAS (N,N) 


GK(N,M) 


GKS(N,M,NSAM) 
H(M,N) 
HS(M,N) 

ITER 

K 


PHI(N,N) 


PHIS(N,N) 


PKK(N,N) 


PKKMi(N,N) 


PKKS(N,N) 


Q(N,N) 


Sample prediction miss distance vector 


covariance of forcing inputs matrix (double 
precision) 


identity matrix (double precision) 
sample estimation error vector 


the sample filter gain matrix (double 
precision) 


the matrix which relates the forcing inputs 
to the state vector (double precision) 


the same as GAMMA but single precision 


the sample filter gain matrix (single 
precision) 


the fiiter gain schedule 

the observation matrix (double precision) 
the observation matrix (single precision) 
the ensemble member index 

the sample index 


the state transition matrix (double 
precision) 


the state transition matrix (single 
precision) 


tne ^“сосуатгтапсе Чо ио с ЕИО ЧЕСТО maLrıxz 
(double precision) 


the one step prediction error covariance 
matrix (double precision) 


iheoewovarqsanececewlwWwadstinsemen error matrix 
(single precision) 


the state excitation covariance matrix 
(double-precision) 
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R(N,N) 


SIGV(N) 


SIGW(N) 


SIGXZ(N) 
V(N) 

W(N) 

X(N) 
XHKK(N) 
XHKKMI(N) 
XM(N,NSAM) 


XS(N,NSAM) 


XZMEAN(N) 


Z(N) 


the covi;.a.--^of measurement noise matrix 
(double ¡res c. ion) 


the measurement noise standard deviation 
vector 


the random tcrcing input standard deviation 
vector 


the initial semple point standard deviation 
the measurement noise vector 

the forcing input vector 

ПО саке CO DS 

the state asstimate vector 

the one step state prediction vector 

the mean st2fo.trajectory over the ensemble 


the state trajectory for the current member 
of the ensemble 


wae mean ОЕ the anitial state samole 


the measurement vector 


A.5 A List of IBM System/360 Source Library Subroutines 


Used by the Monte-Carlo Program 


GAUSS3 
OVFLOW 
PLOTP 


SNORM 


matrix inversion by the method of Gauss 
utility routine required by SNORM 
plot on printer-multiple curves 


normally distributed pseudo-random 
deviations (shuffled) 
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APPENDIX B 


PREDICTION ACCURACY CALCULATION 


The accuracy of the predictions which can be made from 
each sample estimate are of prime concern in a gun fire con- 
trol system since the gun is aimed at a predicted target 
position. A relatively simple procedure for measuring the 
prediction accuracy over realistic time intervals is de- 
scribed below. 

The time interval over which the predictions are made is 
based on a straight line approximation of the time of flight 
(TOF) versus range curve for a 5"/54 caliber projectile [6] 
as indicated in Fig. 4: 

Measurement of the slope oí the straight line approxima- 
tion in Fig. 4 shows that 


TOF 4g 49107 (s/yd)x Range (yds) (1) 


Taking the maximum effective range of the gun against an 
air target into consideration, the time interval over which 
predictions are made is limited to a maximum of ten seconds. 
This corresponds to a maximum effective range of approxi- 
mately 7,000 yards. 

A further approximation iS made in order to simplify the 
matrix manipulations required. The TOF is truncated to an 
integer multiple of the sample period. The prediction is 
based on the equation 


x(k/k-N) = 9$ (Oc-N/k-N) (2) 
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where N is the integer num > ©; sanple periods over which 


» -— 4 
„n 
~ а = 


the prediction is made. 299 
The prediction miss distance is calculated by finding 
the square root of the sum di e >quares of the difference 


between the predicted position states and the respective 


true position states. 


i.e. 


(miss distance)? [X; (k/k-N)-x, (k)]? 


+ 


[X, (k/k-N)-x, (k)]? 


+ 


[£, (k/k-N)-x, (k)]? (3) 


The subroutine which calenlates the prediction miss dis- 
tance in the Monte-Carlo К с stores a running sum of 
miss distances for each sample point over the ensemble so 
that the mean miss distance can be calculated. The subrou- 
tine also counts the number of predictions made (the number 
of samples which meet the projectile firing criteria), the 
number of predictions with miss distance less than 24 feet 
and the number of predictions with miss distances less than 
48 feet for each run over the ensemble. The number of pre- 
dictions with miss distance less than 48 feet at each sample 
point over the ensemble is also recorded. 

ine Fortran listing of subroutine Burst We included as 
part of the Monte-Carlo program in Appendix A. Notice that 
equation (2) is not evaluated in the Computer by raising the 
na matrix Ф to the Nie power by direct matrix algebra, 


which would be very time consuming, but is calculated by 


TOO 





Operating on the individual components of 6. The N-step 


state transition matrix is given by 


om © 0 
N _ 
DEL ONE. 9 
9 DN 


where for the constant-velocity model 


1.0 NT 


0.0 1598 0, 


and for the correlated random acceleration model 


N .n-1 N-1 n=) 
0 NT en 933 +T0,3,2, [N-n]033 
N nen 
9, = | 0.0 10 0,3, Í, 033 


where 
О зч 2 ғ + er am 
POE [1-e 2 
644 = E 


> 


(5) 


(6) 


(Т) 
(8) 


(9) 


When the prediction is to be based on an assumption of 


constant velocity and the estimator is based on the model 


N 
which assumes correlated random accelerations, 0 


ms 


is calcu- 


lated from equation (6) by setting %,3, %23, and %33 equal 


Lo zero. 


Equations (5) and (6) are derived by self-multiplying 


Lhe matrices of equations (3.1) and (3.17) respectively and 


establishing a pattern. 
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APPENDIX C 


A TARGET TRAJECTORY GENERATING PROGRAM 


eser i ptr rono the Program 


In order to test the filter's performance against maneu- 
vering targets it is desirable to have a realistic trajectory 
for which position, velocity and acceleration states are 
known at each sample point. 

The maneuver decided upon is a periodic, planar maneu- 


ver described by the equations 


y(x,t) = pAsin’(kx), (1) 

where 
T O KX S 

and 

e 2 € A 

Piene distance traveled along the EP 
uring the maneuver 

x x E-—IPDSPALXCOSMHXxX)OSihCkX) (2) 

апа 


y(x,X,X,t) = 2pAklxcos(kx)sin(kx) 
+ kx*[2cos*(kx)-1]) (3) 
The "form of the maneuver is shown in Fig. 5. The variables 
765 х, апа х have been used to represent x(t), er) апа x(t) 
“respectively in equations (1) through (14) to make the nota- 


biome mess awkward. 
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The two dimensional geometry of the maneuver yields the 


relationship 
Ve 877 (4) 
where v, the speed of the aircraft tangent to the x-y plane, 
is assumed to be constant. 
Transposing the terms of equation (4), substituting the 
right side of equation (2) for y. and taking the square root 
of both sides yields 


° V 
ll I EA —- —- = O—s А. 
[1*4A?k?cos^(kx)sin^(kx)]? ' | (5) 


Equation (5) is differentiated with respect to time with the 


result 


ү 4A?k?vxsin(kx)cos(kx)[1-2cos?(kx)] 


= 7 (6) 
(ІРА k cos (kx)sin (xx WE 
If the system states are defined as 
Xa. X(t) (7) 
x2 = x(t) (8) 
хз = у(х,%) (9) 
ба - 7050) (10) 
the resulting nonlinear state equations are 
Xi 7 X2 (11) 
оты 4A?k?vx5sin(kxi)cos(kxi) [1-2cos^(kx1)] (12) 
[1*4A?k?cos? (kx, )sin?(kx, ) P/? 
D NT x. (13) 


Xy = 2pAk([x,cos(kx, )sin(kx, )-kx2 ([1-2cos? (kx, ) ]) (14) 


TIS 








The state equations are used in a computer program to 
make up the maneuvering portion of the trajectory. The com- 
puter program utilizes subroutine RKLDEQ, a Runge-Kutta-Gill 
differential equation solver, which is a member of the IBM 
system 360 Source Library. 

In order that filter's performance against a target 
which maneuvers after the filter gains have reached steady 
state can be studied, provision is made for a constant-veloc- 
ity section of target trajectory preceding the maneuver. A 
constant-velocity section of track following the maneuvers 
allows the filter's recovery time to be examined. 

Although the trajectory is originally generated with 
the x-axis as its axis of symmetry, it may be rotated to any 
bearing angie and any iirst quadrant elevation angle. The 
geometry involved in the rotation of the axis of symmetry is 
i hic a TCO Fig. 6 through Fie” 8. As indicated an Laca 
6 for the elevation angle rotation, a new coordinate system 
is established which places the axis of symmetry, now labeled 
the x'-axis, at the desired elevation angle, E. The x and 


z-coordinate states can be calculated from the x' coordinate 


states. 
x(k) = x'(k)cos(E) (15) 
x(k) = x'(k)cos(E) (16) 
X(k) 7» X'(k)cos(E) (17) 
z(k) = x'(k)sin(E) (18) 
z(k) = x'(k)sin(E) (19) 
Z(k) = X'(k)sin(E) (20) 





X(k)2X'(k)cos(E) 


! | X'(k) z(k)=z'(k)sin(E) 
Z'(k)s0 


Axis of 
Symmetry 
2(k)=x'(k)sin(E) 


—— (EEE (dm č — —— —Á— A r e n E 


E=Angle of 
Rotation in 
Elevation 


-e o C s ee «нысы» ë 


X(k)=x'(k)cos(E) 


FIGURE 6 ROTATION OF THE TRACK IN ELEVATION (Vertical 
Plane) 
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FIGURE 7 ROTATION OF THE TRACK AXIS IN BEARING 
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The y-coordinate states are not affected by the elevation 
angle rotation since the rotation is about the y-axis. The 
bearing angle rotation is depicted in Fig. 7 -- the case 
shown is for a zero elevation angle. The axis of symmetry 
is rotated through the angle B, which results in a point on 
the trajectory curve being displaced from the x-axis by a 
total angle of B+ 0. Figure 7 shows the relationship be- 
tween the x- and y-direction velocities before and after 
rotation. The z-coordinate states are not affected by the 
bearing rotation since the rotation is about the z-axis. 
The values of the x and y-coordinate states after rotation 
of the trajectory's axis of symmetry through elevation 


angle E and bearing angle B are given by 


x(k) = R(k)cos[B+9(k)]cos(E) , (53) 

x(k) = [x'(k)cos(B)-y'(k)sin(B)]cos(E) , (22) 

x(k) = [x'(k)cos(B)-y'(k)sin(B)]cos(E) , ce 

y(k) = R(k)sin[B+6(k)]cos(E) , (24) 

y(k) = &'(k)cosEsinBt¥'(k)cosB , 025) 

y(k) = x'(k)cosEsinBty'cosB , (26) 
where 

R(k) = Vx'(k)?+y'(k)? = /x(k)? *y(K)?*z(k)? (27) 
and 

6(k) = arctan [y'(k)/x'(k)] (28) 


Equations (22) and (25) are derived from the trigono- 
m@trice relationship»shown in Fig. Y. Equations (23) and (26) 
are derived from a similar relationship between the accel- 
eration states. 
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Provision is made for the trajectory to have non-zero 
intercepts with the y and/or z axis. This is accomplished 
simply by adding any desired value to the y-and/or z-position 
States at each sample point. 

A subroutine is provided which performs conversion of 
all state values from cartesian to spherical coordinate 
quantities for cases in which spherical coordinate data is 
required. The spherical coordinate data is in units of 


yards in range and radians in bearing and elevation. 


C.2 Relationship Between Velocity, Period of Maneuver, 
Amplitude of tianeuver, and Acceleration 
The magnitude of the total acceleration for the x-y 


plane maneuver is given by 

= /&? 4g? (29) 
where x and y are given by equations (6) and (3) respective- 
ly. Making the substitutions as indicated in equation (29) 
and laboriously solving for a yields 


LANZ COS M (30) 
1*t4A^k?cos?(kx)sin?(kx) 


which can be further simplified by utilizing the trigono- 


Шеше аш оп&м1р 


cos2a - 2cos^a-1 (31) 
and 


ѕіп2а = 2sinacosa (92) 
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with the result 


_ _~2Ak* v*cos(2kx) 


(33) 
1 + A*k*sin’ (2kx) 


To find the point, with respect to x-position, at which the 
maximum acceleration occurs during the maneuver it is nec- 
essary to differentiate equation (33) with respect to x and 


set the result equal to zero. This yields 
4A*k* v7 sin(2kx)[1+4Ak* sin” (2kx)+Ak* cos? (2kx)] = 0 (34) 
which is true when 


sin(2kx) - 0, 


ЛЕР 2Ех-П,) п=0,‚1,‚2,.... 
Substitution of 2kx - nm into equation (28) yields 


[а | = 2Ak^v^* (35) 
max 


Equation (35) gives the relationship between target speed, 
maneuver amplitude, maneuver period, and the maximum accel- 
eration required to complete the maneuver. This relation- 
ship Is useful as an aid in choosing values for the variables 


involved. 


O ample Saci wl mack. Generat 1nggProgsram Use 
Example LEA trajectory 15 £o be generated which ap- 


Proacheswriiesporıecın Troms aerange201 3500 y2rQ4s!with au 20 
degree dive angle at a constant speed of 200 yd/s. А 2.5 g 
maneuver of approximately ten seconds duration is to com- 
mence when the target is at a range of 2,700 yards. Sixty 


samples are to be generated at a four-hortz sampling rate. 
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Телгә task is) te ne track specifücations 


into theBdaralrequired Ly e generating program. 


A maneuver duration of ten seconds is specified. To 


calculate the distance таушы ыйїопк the axis of symmetry 


during a 2.5 g, ten second; sin? maneuver would not be a 


simple task. However expe- aws i has shown that the change 
in range rate during such a maneuver is relatively smalls, 


therefore the distance is approximately 


х cr = (10s)(200yd/s) = 2000 yd. (37) 


Equation (36) can now be used to solve for the amplitude of 


the maneuver. 


SS 


-_2.5(g) x 10.73(yd/s.g) 


2 oa) 1200(ya/s)7? 


(38) 


35107 Vd. 


The four hertz sampling rate inverts to a 0.25 s sampling 
period. 

The data is now available in the proper form to make up 
the data input cards for the track generating program. The 


cards must contain the following information. 


Card 1: 
Б СЕ Initia lerange = 3500.0 
S(2,1) = initial range rate = -200.0 


Toi 








Card 2: 


RATE = sampling period = 0025 
V — target speed = 200.0 
A = maneuver amplitude = 33.290 
XPER = maneuver duration = 2000.0 
BR = bearing angle = 0.0 (no bear- 
ing specified) 
E = elevation angle = 20.0 
RMAN = maneuver commencement = 2700.0 
range 
TRKEND = trajectory end range =-9999.9 
ë (the track is to be ended 
by the sample count) 
Band 3: 
CPA = y-direction displacement = 0.0 
HO = z-direction displacement = B 
Card 4: 
NSAM = no. of samples desired = 60 
N — no. of states generated = 9 
ICOORD = coordinate system indi- = 0 
cator (data 1s tompemin 
cartesian coordinates) 
ПРОМОВ т саса output indicator - i 


(the output data is to 
be punched onto cards) 


Table VIII lists the data generated by the track gener- 
ating program for this example. 

EAU a G, e R a a a e 
ginning at the point x=5000 yards, y=1000 yards and z=333.3 
yards in cartesian coordinates and proceeding in the negative 
x-direction with a speed of 200 yds/s until the x-direction 
displacement is 2800 yards is to be produced. The sampling 
period is 0.25 seconds. The output data is to be in spheri- 
cal coordinates. 

This trajectory can be generated quite simply by uti- 
lizing the CPA and HO features of the program. Since tho 


trajectory is generated initially with the x-axis as its 
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TRACK GENERATING PROGRAM OUTPUT FOR EXAMPLE 1 


TABLE VIII 
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axis of symmetry and the y and z-direction displacements 

are added afterward, the specified values of initial x-posi- 
tion and velocity are taken as the initial range and range- 
Pre respectively. The other required input values are 
given directly in the track specifications. 


The input data cards contain the following information. 


Card 1: 
S(1,1) = initial range (x-position) = 5000.0 
S(2,1) = initial range rate (x-velocity) = -200.0 
care: 
RATE = sampling period = 0.25 
V = target speed = 200.0 
A = maneuver amplitude = 0.20 
XPER = maneuver duration = 0.0 
BR = bearing angle = 0.0 
E = elevation angle = 0.0 
RMAN = maneuver commencement range = 0.0 
TRKEND = trajectory end range = 2800.0 
Cara 3E 
CPA — y-direction displacement = 1000.0 
HO = z-direction displacement "2999 
Card 4: 
NSAM = number of samples = 200 
(the track is range limited) 
N = number of states generated = 9 
ICOORD = coordinate system indicator = 1 
(spherical coordinate data desired) 
TPBUNCAS="ea rd output indicator = О 


(по cards are desired) 
Table IX lists the trajectory data generated by the 


track generating program for this example. 
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TABLE IX TRACK GENERATING PROGRAM OUTPUT FOR EXAMPLE 2 


OUTPUT DATA IN SPHERICAL COORDINATES 


R2DOT BEARING BDOT B2COT ELEVATION EDOT E2DCT 
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APPENDIX D 


ADAPTIVE FILTER SUBROUTINES 


FORTRAN listings of the three adaptive filter subrou- 
tines are included in this appendix. Also included is a 
flowchart which is based on the first version of the switch- 
on-residuals adaptive filter subroutine, because there is 
essentially no difference in the program flow, there is no 


separate flowchart provided for the second version. 
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